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Abstract 

In this work we introduce the Multi-Index Stochastic Collocation method (MISC) for computing statistics 
of the solution of a PDE with random data. MISC is a combination technique based on mixed differences 
of spatial approximations and quadratures over the space of random data. We propose an optimization 
procedure to select the most effective mixed differences to include in the MISC estimator: such optimization 
is a crucial step and allows us to build a method that, provided with sufficient solution regularity, is 
potentially more effective than other multi-level collocation methods already available in literature. We 
then provide a complexity analysis that assumes decay rates of product type for such mixed differences, 
showing that in the optimal case the convergence rate of MISC is only dictated by the convergence of the 
deterministic solver applied to a one dimensional problem. We show the effectiveness of MISC with some 
computational tests, comparing it with other related methods available in the literature, such as the Multi- 
Index and Multilevel Monte Carlo, Multilevel Stochastic Collocation, Quasi Optimal Stochastic Collocation 
and Sparse Composite Collocation methods. 

Keywords: Uncertainty Quantification, Random PDEs, Multivariate approximation, Sparse grids, 
Stochastic Collocation methods, Multilevel methods, Combination technique. 

2010 MSC: 41A10, 65C20, 65N30, 65N05 


1. Introduction 

Uncertainty Quantification (UQ) is an interdisciplinary, fast-growing research area that focuses on de¬ 
vising mathematical techniques to tackle problems in engineering and natural sciences in which only a 
probabilistic description of the parameters of the governing equations is available, due to measurement er¬ 
rors, intrinsic non-measurability/non-predictability, or incomplete knowledge of the system of interest. In 
this context, “parameters” is a term used in broad sense to refer to constitutive laws, forcing terms, domain 
shapes, boundary and initial conditions, etc. 

UQ methods can be divided into deterministic and randomized methods. While randomized techniques, 
which include the Monte Carlo sampling method, are essentially based on random sampling and ensemble 
averaging, deterministic methods proceed by building a surrogate of the system’s response function over the 
parameter space, which is then processed to obtain the desired information. Typical goals include computing 
statistical moments (expected value, variance, higher moments, correlations) of some quantity of interest of 
the system at hand, typically functionals of the state variables (forward problem), or updating the statistical 
description of the random parameters given some observations of the system at hand (inverse problem). In 
any case, multiple resolutions of the governing equations are needed to explore the dependence of the state 
variables on the random parameters. The computational method used should therefore be carefully designed 
to minimize the computational effort. 
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In this work, we focus on the case of PDEs with random data, for which both deterministic and ran¬ 
domized approaches have been extensively explored in recent years. As for the deterministic methods, 
we mention here the methods based on polynomial expansions computed either by global Galerkin-type 
projections [1, 2, 3, 4, 5] or collocation strategies based on sparse grids (see e.g. [6, 7, 8, 9]), low-rank 
techniques [10, 11, 12, 13] and reduced basis methods (see e.g. [14, 15]). All these approaches have been 
found to be particularly effective when applied to problems with a moderate number of random parameters 
(low-dimensional probability space) and smooth response functions. Although significant effort has been 
expended on increasing the efficiency of such deterministic methods with respect to the number of random 
parameters (see, e.g., [16], the seminal work on infinite dimensional polynomial approximation of elliptic 
PDEs with random coefficients), Monte Carlo-type approximations remain the primary choice for problems 
with non-smooth response functions and/or those that depend on a high number of random parameters, 
despite their slow convergence with respect to sample size. 

A very promising methodology that builds on the classical Monte Carlo method and enhances its perfor¬ 
mance is offered by the so-called Multilevel Monte Carlo (MLMC). It was first proposed in [17] for applica¬ 
tions in parametric integration and extended to weak approximation of stochastic differential equations in 
[18], which also provided a full complexity analysis. Let {hg}f =0 be a (scalar) sequence of spatial/temporal 
resolution levels that can be used for the numerical discretization of the PDE at hand and {Ff\^ =0 be the 
corresponding approximations of the quantity of interest, and suppose that the final goal of the UQ analysis 
is to compute the expected value of F, E[F], While a classic Monte Carlo approach simply approximates 
the expected value by using an ensemble average over a sample of independent replicas of the random 
parameters, the MLMC method relies on the simple observation that, by linearity of expectation, 

L 

E[F] » E[F l ] = E[F 0 ] + ^E[F* - -Flg-i], (1) 

t =i 

and computes by independent Monte Carlo samplers each expectation in the sum. Indeed, if the discretiza¬ 
tion of the underlying differential model is converging with respect to the discretization level, £, the variance 
of (Ft — Fe_i) will be smaller and smaller as £ increases, i.e., when the spatial/temporal resolution increases. 
Dramatic computational saving can thus be obtained by approximating the quantities E [Fg — Fg_ i] with a 
smaller and smaller sample size, since most of the variability of F will be captured with coarse simulations 
and only a few resolutions over the finest discretization levels will be performed. The MLMC estimator is 
therefore given by 


L 


E[F]«£ 

£=0 


1 M e 

(Fg(u! m ^) — Fg-i(co m ,e)) > 

^ m= 1 


with F_ i(-) = 0, 


( 2 ) 


where are the i.i.d. replicas of the random parameters. The application of MLMC methods to UQ 
problems involving PDEs with random data has been investigated from the mathematical point of view in 
a number of recent publications, see e.g. [19, 20, 21, 22, 23]. Recent works [24, 25, 26, 27] have explored the 
possibility of replacing the Monte Carlo sampler on each level by other quadrature formulas such as sparse 
grids or quasi-Monte Carlo quadrature, obtaining the so-called Multilevel Stochastic Collocation (MLSC) or 
Multilevel Quasi-Monte Carlo (MLQCM) methods. See also [28] for a related approach where the Multilevel 
Monte Carlo method is combined with a control variate technique. 

The starting point of this work is instead the so-called Multi-Index Monte Carlo method (MIMC), recently 
introduced in [29] , that differs from the Multilevel Monte Carlo method in that the telescoping idea presented 
in equations (l)-(2) is applied to discretizations indexed by a multi-index rather than a scalar index, thus 
allowing each discretization parameter to vary independently of the others. Analogously to what done in 
[24, 25, 26] in the context of stochastic collocation, here we propose to replace the Monte Carlo quadrature 
with a sparse grid quadrature at each telescopic level, obtaining in our case the Multi-Index Stochastic 
Collocation method (MISC). In other words, MISC can be seen as a multi-index version of MLSC, or a 
stochastic collocation version of MIMC. From a slightly different perspective, MISC is also closely related to 
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the combination technique developed for the solution of (deterministic) PDEs in [30, 7, 31, 32, 33]; in this 
work, the combination technique is used with respect to both the deterministic and stochastic variables. 

One key difference between the present work and [24, 25, 26] is that the number of problem solves to be 
performed at each discretization level is not determined by balancing the spatial and stochastic components 
of the error (based, e.g., on convergence error estimates), but rather suitably extending the knapsack- 
problem approach that we employed in [34, 35, 36] to derive the so-called Quasi-Optimal Sparse Grids 
method (see also [37]). A somewhat analogous approach was proposed in [38], where the number of solves 
per discretization level is prescribed a-priori based on a standard sparsification procedure (we will give more 
details on the comparison between these different methods later on). In this work, we provide a complexity 
analysis of MISC and illustrate its performance improvements, comparing it to other methods by means of 
numerical examples. 

The remainder of this paper is organized as follows. In Section 2, we introduce the problem to be 
solved and the approximation schemes that will be used. The Multi-Index Stochastic Collocation method 
is introduced in Section 3, and our main theorem detailing the complexity of MISC for a particular choice 
of an index set is presented in Section 4. Finally, Section 5 presents some numerical tests, while Section 6 
offers some conclusions and final remarks. The Appendix contains the technical proof of the main theorem. 
Throughout the rest of this work we use the following notation: 

• N denotes the set of integer numbers including zero; 

• N + denotes the set of positive integer numbers, i.e. excluding zero; 

• R + denotes the set of positive real numbers, R + = {r G R : r > 0}; 

• 1 denotes a vector whose components are always equal to one; 

• ef denotes the £-th canonical vector in R K , i.e., (ef)i = 1 if i = i and zero otherwise; however, for 
the sake of clarity, we often omit the superscript k when obvious from the context. For instance, if 
v G R w , we will write v — e.\ instead of v — e^; 

• given v £ R N , |w| = J2n =l v m max ( u ) = max„ = i i ...jv v n and min(u) = min„ =lj ..,jv v n ; 

• given v £ M. N and / : R —► R, f{v) denotes the vector obtained by applying / to each component of 

V, f(v) = [f(v- l), f(v 2 ), • • • , f(v at)] G R^; 

• given v, w £ R N , the inequality v > w holds true if and only if v n > w n Vn = 1,..., N. 

• given v £ R D and w £ R w , [v, w] = (iq,..., vd,w±, ..., wn) G R D+Ar . 

2. Problem setting 

Let 23 C R d , d, = 1,2,3, be an open hyper-rectangular domain (referred to hereafter as the “physical 
domain”) and let y = (j/i, 2 / 2 5 ■ • ■,2/jv) be a TV-dimensional random vector whose components are mutually 
independent and uniformly distributed random variables with support T n C R and probability density 
function p n {yn ) = . Denoting r = Tixr 2 ...xTjv (referred to hereafter as the “stochastic domain” 

or “parameter space”) and by ctb(T) the Borel cr-algebra over T, p(y)dy = n^=i Pn{yn)dy n is therefore a 
probability measure on T, due to the independence of y n , and (T, ctb(T), p(y)dy) is a complete probability 
space. Consider the following generic PDE, together with the assumption stated next: 

Problem 1. Find u : 23 x T —> R such that for p-almost every y G T 

jC(u; x, y) = T(a;) x G 23, 

|it(a;, y) = 0 x £ 923. 

Assumption 1 (Well posedness). Problem 1 is well posed in some Hilbert space V for p-almost every y G T. 
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The solution of Problem 1 can be seen as an iV-variate Hilbert-space valued function u{y) : P —> V. 
The random variables, y n , can represent scalar values whose exact value is unknown, or they can stem 
from a spectral decomposition of a random field, like a Karhunen-Loeve or Fourier expansion, possibly 
truncated after a finite number of terms, see, e.g., [6, 36]. It is also useful to introduce the Bochner space 
L 2 p (T; V) = {u : T —>■ V strongly measurable s.t. f r \\u(y)\\y p(y)dy < oo} . Finally, given some functional 
of the solution u, 0 : V —> K, we denote by F : T —> K. the iV-variate real-valued function assigning to each 
realization y G F the corresponding value of 0[it] (quantity of interest), i.e., F(y ) = 0[u(-,y)], and we aim 
at estimating its expected value, 

E [F) = j F{y)p(y)dy. 

Example 1. As a motivating example, consider the following elliptic problem: find u : 23 x F —> K. such 
that for p-almost every y GT 

f-di v(a(x,y)Vu(x,y)) =3 r (x) x G 

y) = h(x) x G 923, 

holds, where div and V denote differentiation with respect to the physical variables, x, only, and the function 
a : 23 x F —► R. is bounded away from 0 and oo, i.e., there exist two constants, a m i n , a ma x! such that 

0 < a m i n < a(x, y) < a max < oo, \/x G 23 and for p-almost every y GT. (4) 

This boundedness condition guarantees that Assumption 1 is satisfied, i.e. the equation is well posed for 
p-almost every y G T, thanks to a straightforward application of the Lax-Milgram lemma; moreover, the 
equation is well posed in L 2 (T;V), where V is the classical Sobolev space Ft q(23), see, e.g., [6]. This is the 
example we will focus on in Section 5, where we will test numerically the performance of the Multi-Index 
Stochastic Collocation method that we will detail in Section 3. 

Remark 1. The method that we present in the following sections can be also applied to more general 
problems than Problem 1 in which the forcing terms, boundary conditions and possibly domain shape are 
also modeled as uncertain; the extension to time-dependent problems with uncertain initial conditions is 
also straightforward. Other probability measures can also be considered; the very relevant case in which the 
random variables, y n , are normally distributed is an example. 

Remark 2. As will be clearer in a moment, the methodology we propose uses tensorized solvers for deter¬ 
ministic PDEs. Although for ease of exposition we have assumed that the spatial domain, 23, is a hyper¬ 
rectangle, it is important to remark that the methodology proposed in this work can also be applied to non 
hyper-rectangular domains: this can be achieved by introducing a mapping from a reference hyper-rectangle 
to the generic domain of interest (with techniques such as those proposed in the context of Isogeometric 
Analysis [39] or Transfinite Interpolation [40]) or by a Domain Decomposition approach [41] if the domain 
can be obtained as a union of hyper-rectangles. 

2.1. Approximation along the deterministic and stochastic dimensions 

In practice, we can only access the value of F via a numerical solver yielding a numerical approximation 
of the solution u of Problem 1, which depends on a set of D discretization parameters, such as the mesh-size, 
the time-step, the tolerances of the numerical solvers, and others, which we denote by hi,i = 1,... ,D\ we 
remark that in general D , the number of parameters, might be different from d, the number of spatial 
dimensions. For each of those parameters, we introduce a sequence of discretization levels, h i a , a = 1, 2,..., 
and for each multi-index a G N5, we denote by u a (x,y) the approximation of u obtained from setting 
hi = hi t0li , with the implicit assumption that u a {x,y) —> u(x,y) as mini <i<DOn —> oo for p-almost every 
y G T; similarly, we also write F a (y) = 0[u“(-, y)\. For instance, we could solve the problem stated in 
Example 1 by a finite differences scheme with grid-sizes = ho2~ ai in direction i = 1,..., D, for some 
ho > 0. 
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The discretization of F a over the random parameter space T will consist of a suitable linear combination 
of tensor interpolants over T based on Lagrangian polynomials. Observe that this approach is sound only 
if F a is at least a continuous function over T (the smoother F a is, the more effective the Lagrangian 
approximation will be); for instance, for the problem stated in Example 1, it can be shown under moderate 
assumptions on a(x,y) that F and F a are y-analytic, see, e.g., [35, 16]; we will return to this point in 
Section 5. 

To derive a generic tensor Lagrangian interpolation of F a , we first introduce the set C°(r n ) of real-valued 
continuous functions over T n , and the subspace of polynomials of degree at most q over L n . P 9 (T ra ) C C°(T rl ). 
Next, we consider a sequence of univariate Lagrangian interpolant operators in each dimension Y n , i.e., 
{U™^"^}/3„eN + ; where we refer to the value j3 n as the “interpolation level”. Each interpolant is built over 
a set of m(/3 n ) collocation points, = {y^, } C T n , where to is a strictly increasing 

function, with m(0) = 0 and to( 1) = 1, that we call the “level-to-nodes function”; thus, the interpolant 
yields a polynomial approximation, 

m(/3„) / m(t 3„) _ k 

U™^[f](y n ) = Y, \ fivi) II 

j=i V fc=i,**y Vn ~ y » 

with the convention that [/] = 0 V/ £ C°(r„). 

The 7V-variate Lagrangian interpolant can then be built by a tensorization of univariate interpolants: de¬ 
note by C°(T) the space of real-valued N-variate continuous functions over T and by P' ? (T) = ® n=1 P 9 "(r n ) 
the subspace of polynomials of degree at most q n over T n , with q = (qi, ..., qjv) £ N^, and consider a 
multi-index /3 € N* assigning the interpolation level in each direction, y n \ the multivariate interpolant can 
then be written as 

U m(/3) . C 0 (p) pM-i( I ') ] U m(/3) [E“](y) = (ll™ (/3l) 0 • • • ® [F a ](y). 

The set of collocation points needed to build the tensor interpolant IP"^ [it] (y) is the tensor grid = 

with cardinality #T m ^) = nti m(P n ). Observe that the Lagrangian interpolant immediately 
induces an iV-variate quadrature formula, : C°(T) —> R, 



Q m(/3)[ F «] =E \u m ^[F a ](y) 


|j"(W 

Y pa (y3) w n 


3=1 


where yj £ and the quadrature weights vjj are the expected values of the Lagrangian polynomials 

centered in yj , which can be computed exactly for most of the common interpolation knots and probability 
measures of the random variables. 

It is recommended that the collocation points J~C n ' to be used in each direction are chosen according 
to the underlying probability measure, p(y n )dy n , to ensure good approximation properties of the interpolant 
and quadrature operators, IP "^ and Q m ^P\ Common choices are Gaussian quadrature points like Gauss- 
Legendre for uniform measures or Gauss-Hermite for Gaussian measures, cf. e.g., [42], which are however 
not nested, i.e., This means that they are not optimal for successive refinements of 

the interpolation/quadrature, and we will not consider them in this work. Instead, we will work with nested 
collocation points, and specifically with Clenshaw-Curtis points [34, 43], that are a classical choice for the 
uniform measure that we are considering here; other choices of nested points are available for uniform random 
variables, e.g., the Leja points [34, 44], whose performance is somehow equivalent to that of Clenshaw-Curtis 
for quadrature purposes, see [45, 46]. Clenshaw-Curtis points are defined as 

yi=cosf ^ 1 ' >7r \ , 1 <j<m(i n ), (5) 

\m(i n ) -1 J 

together with the following level-to-nodes relation, m(i n ), that ensures their nestedness: 

m(0) = 0, m(l) = 1, m(i n ) = 2 1 " -1 + 1. 
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( 6 ) 








We conclude this section by introducing the following operator norm, which acts as a “Lebesgue constant' 
from C°(r) to L 2 p {Y)\ 


N 


•d/3) 


= J[[M^" ) , with sup WK^fh^y 


lli°°(r„) = 


r f~Y~L f [3 1 

In particular, for the Clenshaw-Curtis points, it is possible to bound M n K J as: 


N 


ir(/3) < m ^(/ 3) = p 


ff™(/3n) 

11 n.est i 


M* t = 


n—1 


for 9 = 1 


— log(g — 1) + 1 for q > 2. 

7T 


(7) 


( 8 ) 


See [34] and references therein. 

Remark 3. Nested collocation points have been studied also for other probability measures than uniform 
probability measures. In the very relevant case of a normal distribution, one possible choice is the Genz- 
Keister points [47, 36]; we mention also the recent work [46] on generalized Leja points that can be used for 
arbitrary measures on unbounded domains. 

3. Multi-Index Stochastic Collocation 

It is easy to see that an accurate approximation of E [F] by a direct tensor technique as the one just 
introduced, E[F] ss Q m (/3)[jr a ] ; might require a prohibitively large computational effort even for moderate 
values of D and N (what is referred to as the “curse of dimensionality”). In this work, following the setting 
that was presented in [34, 29], we propose the Multi-Index Stochastic Collocation as an alternative. It can 
be seen as a generalization of the telescoping sum presented in the introduction, see equations (1) and (2). 
Denoting Q m (P )(F“] = the building blocks of such a telescoping sum are the first-order difference 

operators for the deterministic and stochastic discretization parameters, denoted respectively by A(* et with 
1 < * < D and Af oc with 1 < j < N: 


A t et \F a .p\ 

&f oc [F a , 0 \ 


Fct,(3 Fat—eii/3-) 

if ai > 1, 

(9) 

{Fol,{3 

if ai = 1, 

\ ^cx,f3 Fot,(3—ej i 

if dj > 1, 

(10) 


if fy = 1. 


We then define the first-order tensor difference operators, 


D 

A det [F a ,p] = 0 A i et [F a , p \ = A det [ A det [ • • • A det [F a , p ]) ] 

N 

A st ° c [F a ,p} = 0 Af-[F„ i/3 ] = Y, ("I ■ 
j = 1 je{o,i}« 


£ (-1 )WF a - jifi , 

je{ o,i} D 


( 11 ) 

( 12 ) 


with the convention that F v w = 0 whenever a component of v or w is zero. Observe that computing 
A det [F ajl g] actually requires up to 2° solver calls, and analogously applying A stoc [f 1 aj| g] requires interpo¬ 
lating F a on up to 2 n tensor grids; for instance, if D = N = 2 and a, [3 > 1, we have 

A det [ F a , p ] = A det [ A det [ F a p} ] = A det {Fol.p - F a _ eit/3 ] = F ai/3 - F a _ eii(9 - F Q _ e2i(9 + F a _ ltf3 , 

A & t°c[F a ^g] = F a ^p — -F a ,/3- ei — F a {3— e2 + F a ^ 3_i. 
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Finally, letting A [F a ^\ = A stoc [A dot [F ct;| a]], we define the Multi-Index Stochastic Collocation (MISC) 
estimator of E [F] as 

M : x[F}= Y A 1 F <*a) = E (13) 

[a,/3]e I [a,/3]el 

where X C N^ +Ar and c Qj /3 £ Z. Observe that many of the coefficients in (13), c a ,p, may be zero: in 
particular, c a ,p is zero whenever [a,/3] + j £ X Vj £ {0, l} Ar+D . Similarly to the analogous sparse grid 
construction [34, 48, 7], we shall require that the multi-index set X be downward closed, i.e., 


w r ^. 1 ^- 7 - I a - ei £ X for 1 < i < D and cq > 1, 

V \OL. p G -L , \ 

|/3 — ej £ J for 1 < j < TV and /3j > 1. 

Remark 4. In theory, a MISC approach could also be developed to approximate the entire solution u(x, y) 
and not just the expectation of functionals, considering differences between consecutive interpolant operators, 
U m (P'), on the stochastic domain rather than differences of the quadrature operators, Q m ^\ as a building 
block for the A stoc operators, as well as considering the discretized solution u a rather than just the quantity 
of interest, F a , in the construction of the A det operators. 


3.1. A knapsack-like construction of the setX 

The efficiency of the MISC method in equation (13) will heavily depend on the specific choice of the 
index set, X; in the following, we will first propose a general strategy to derive quasi-optimal sets and then 
prove in Section 4 a convergence result for such sets under some reasonable assumptions. 

To derive an efficient set, X, we recast the problem of its construction as an optimization problem, in 
the same spirit of [29, 34, 35, 7, 37]. We begin by introducing the concepts of “work contribution”, AW a ^p, 
and “error contribution”, A E a p, for each hierarchical surplus operator, A[X’ q , i( 3 ]. The work contribution 
measures the computational cost (measured, e.g., as a function of the total number of degrees of freedom, or 
in terms of computational time) required to add A [F a ,p\ to Mi [F], i.e., to solve the associated deterministic 
problems and to compute the corresponding interpolants over the parameter space, cf. equations ( 11 ) and 
(12); the error contribution measures instead how much the error |E[F] — Mx[F]| would decrease once the 
operator A [F a ,p] has been added to Mx[F]. In formulas, we define 

A W a ,p = Work[Mx U{[cti/ 3 ]} [F]] - Work[Mx[F]] = Work[A[F a ,p]], 


so that 

Work[Mx[F]]= Y AW <x,P, (14) 

K/3]ex 


Observe that this work definition is sharp only if we think of building the MISC estimator with an incremental 
approach, i.e., we assume that adding the multi-index ( a,/3) to the index set X would not reduce the work 
that has to be done to evaluate the MISC estimator on the index set. This implies that one cannot take 
advantage of the fact that some of the coefficients in (13), c a ^p, that are non-zero when considering the set X 
could become zero if the MISC estimator is instead built considering the set XU {[a,/3]}, hence it would be 
possible not to compute the corresponding approximations F a p. This approach is discussed in Section 5.3. 

Similarly, we define 


AE a p 


l%U{[a,/3]}[-F] ^ Mx[F] 


|A[F ai/3 ]|. 


Thus, by construction, the error of the MISC estimator (13) can be bounded as the sum of the error 
contributions not included in the estimator Mx[F], 


Error[M X [F]] = |E[F] - M X [F]| = 


E A 

[a,pm 


< Y i a i f «./3]i = E AE <*#■ ( 15 ) 

\ a ,pm [oc.pm 
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Consequently, a quasi-optimal set X can be computed by solving the following “binary knapsack problem^ 
[49]: 


maximize E 

[a,/3]GN+ +JV 

such that E AW a ,0X a> f3 < Wmax, (16) 

[a,/3]eN° +JV 

Xot,j3 £ { 6 , 1 } , 


and setting I = {[a,/3] £ N+ +n ■ Xa,p = 1}- Observe that such a set is only “quasi” optimal since 
the error decomposition (15) is not an exact representation but rather an upper bound. The optimization 
problem above is well known to be computationally intractable. Still, an approximate greedy solution (which 
coincides with the exact solution under certain hypotheses that will be clearer in a moment) can be found 
if one instead allows the variables x at p to assume fractional values, i.e., it is possible to include fractions 
of multi-indices in X. For this simplified problem, the resulting problem can be solved analytically by the 
so-called Dantzig algorithm [49]: 

1. compute the “profit” of each hierarchical surplus, i.e., the quantity 


Pol,13 ~ 


AW 0il j ; 


2. sort the hierarchical surpluses by decreasing profit; 

3. add the hierarchical surpluses to Mz[-F] according to such order until the constraint on the maximum 
work is fulfilled. 

Note that by construction x a ,p = 1 for all the multi-indices included in the selection except for the last one, 
for which x a .p < 1 might hold; in other words, the last multi-index is the only one that might not be taken 
entirely. However, if this is the case, we assume that we could slightly adjust the value W max , so that all 
x a ,/3 have integer values (see also [7]); observe that this integer solution is also the solution of the original 
binary knapsack problem (16) with the new value of W max in the work constraint. Thus, if the quantities 
A E a j3 and AWq^ were available, the quasi-optimal index set for the MISC estimator could be computed 
as 

X = X(e) = {[a,/3]eN? +JV : > ej , (17) 

for a suitable e > 0. 

Remark 5. The MISC setting could in principle include the Multilevel Stochastic Collocation method pro¬ 
posed in [ 24 ] as a special case, by simply considering a discretization of the spatial domain on regular meshes, 
and letting the diameter of each element (the mesh-size) be the only discretization parameter, i.e., D = 1. 

However, the sparse grids to be used at each level are determined in [24] by computing the minimal number 
of collocation points needed to balance the stochastic and spatial error. This is done by relying on sparse grid 
error estimates; yet, since in general it is not possible to generate a sparse grid with a predefined number 
of points, some rounding strategy to the sparse grid with the nearest cardinality must be devised, which may 
affect the optimality of the multilevel strategy. In the present work , we overcome this issue by relying instead 
on profit estimates to build a set of multi-indices that simultaneously prescribe the spatial discretization and 
the associated tensor grid in the stochastic variables. Furthermore, only standard isotropic Smolyak sparse 
grids are considered in the actual numerical experiments in [ 24 ] (although in principle anisotropic sparse 
grids could be used as well, provided that good convergence estimates for such sparse grids are available), 
while our implementation naturally uses anisotropic stochastic collocation methods at each spatial level. 




The MISC approach also includes as a special case the “Sparse Composite Collocation Method” developed 
in [38], by considering again only one deterministic discretization parameter, i.e., D = 1, and then setting 

Z= j[a,/3]eN^ + *:a + £ft,<it;j, (18) 

with w £ N + . In other words, the approach in [38] is based neither on profit nor on error balancing. 

4. Complexity analysis of the MISC method 

In this section, we assume suitable models for the error and work contributions, A E a p and AW a ,p 
(which are numerically verified in Section 5 for the problem in Example 1) and then state our main con¬ 
vergence theorem for the MISC method built using a particular index set, I*, which can be regarded as an 
approximation of the quasi-optimal set introduced in the previous section. 

Assumption 2. The discretization parameters, hi, for the deterministic solver depend exponentially on the 
discretization level a,, and the number of collocation points over the parameter space grows exponentially 
with the level ft: 

hi,ai — M 1 and A m(ft) A C m ^ U p2j^ z . 

Assumption 3. The error and work contributions, AE a p and AW a p, can be bounded as products of two 
terms, 

AE a ^ < AE^AEp oc and AW a>p < AW* et AW s p toc , 

where AW^ et and A .E det denote the cost and the error contributions due to the deterministic difference 
operator, A det [.F ai| g], and similarly AWp oc and AEp oc denote the cost and the error contribution due to 
the stochastic difference operator, A stoc [i 7 ft | g], cf. equations ( 1 1 ) -( 12 ) . 

Assumption 4. The following bounds hold true for the factors appearing in Assumption 3: 


D 


AW det < 

i = 1 


(19) 

D 



A£ det < C' e d r e r t or [](ft ;ai ) r ~ i , 

2=1 


( 20 ) 

N 

N 


AW| toc < C^f k [] m(/ 3 n ) < O 

IT 9 -. 

( 21 ) 

n—1 

n—1 


A pstoc < ^stoc 

^ error ° 5 


( 22 ) 


for some rates 7 i,r^,Qi > 0 . 

With these assumptions, we are now ready to state our main theorem. The proof is technical and 
we therefore place it in the appendix. The proof is based on summing the error contributions outside a 
particular index set, I*, and the work contributions inside the same index set. This can be seen as a 
weighted cardinality argument in finite dimensions. See also [50, 51] for similar arguments for different 
choices of finite and infinite dimensional index sets. 

Theorem 1 (MISC computational complexity). Under Assumptions 2 to ft the bounds for the factors 
appearing in Assumption 3 can be equivalently rewritten as 

A W a>p < C wolk eZ?=^ ai e 5 
AE a<p < C eiIOT e~ "EiLi r i a i e ~ 9j ex p(^ft-), 
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(23a) 

(23b) 


with 7 i = 7 i log 2 7 77 = 77 log 2 , 6 = log 2 and gi = giCm^ow Define the following set 


D 


N 


X*{L) = l [a, 0 \ e N^ +Ar : + $](<% + < L } L e K+. (24) 


T/ien i/iere exists a constant Qw such that, for any W max satisfying 

W max > C^exp (x), 

and choosing L as 


L = L(W max ) = — ( log 
X 


W„ 


■w 


- (3 - l)log ( - log ^ max 


X 




(25) 


(26) 


with S = 


’ ' ' ' ’ ID+rD 

estimator Mz*(L(Wma. x )) satisfies 


), X = max(H). C = min^ and 3 = #{* = !,...£> : ^ = C}> MJ5C7 


Work[Mx*(i(W max ))] < Wmax, 

Err 01 '[ :M I*(L(W^ max ))] 
w max too W m L (log (W max )) (C+ )(3 } 


= 6 e < 00. 


(27a) 

(27b) 


Remark 6 . The set I* proposed in Theorem 1 can he obtained by assuming that the bounds in equations 
(23a) and (23b) are actually equalities and by using the definition of the quasi-optimal set (18): 


( „ xr e -HfLiriOii e -'E?=iBiexp(SP j ) 1 

- ^ : - e E°iT.a, e i|g| - a j 

= I [ol,0\ G N° +n : ^2(r i +'y i )a i + ^2{5/3 i +gie Sf)i ) < l\, 

l i=l i=l J 


where the last equality holds with L = — loge. 

Remark 7. Refining along the spatial or the stochastic variables has different effects on the error of the 
MISC estimator. Indeed, due to the double exponential e~ ^-o =1 9j eyL vAPj) in (23b), the stochastic contribution 
to the error will quickly fade to zero, which in turn implies that most of the work will be used to reduce the 
deterministic error. This is confirmed by the fact that the error convergence rate in Theorem 1 only depends 
on 7 j and ri, i.e., the cost and error rates of the deterministic solver, respectively. This observation coincides 
with that in [38, page 2299]: “since the stochastic error decreases exponentially, the convergence rate should 
tend towards the algebraic rate of the spatial discretization [...]; see Proposition 3.8 V . Compared with the 
method proposed in [38], MISC takes greater advantage of this fact, since it is based on an optimization 
procedure, cf. equation (17); this performance improvement is well documented by the comparison between 
the two methods shown in the next section. Figure 1 shows the multi-indices included in T according to (24) 
for increasing values of L, for a problem with N = D = 1, 74 = 1, r = 2, and g = 1.5: as expected, the 
shape ofT becomes more and more curved as L grows, due to this lack of balance between the stochastic and 
deterministic directions. 


Remark 8 . Theorem 1 is only valid in case Assumptions 1- 4 cine true. In the next section we motivate these 
assumption for a specific elliptic problem that we use for numerically testing the MISC method. However, we 
stress that deriving bounds on the error and work contributions is problem-dependent and the corresponding 
analysis must be carried out in each case. Moreover, under a different set of assumption the complexity 
theorem would have to be rewritten accordingly. 
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Figure 1: Index sets I(L) for D = N — 1, according to equation (24). 


5. Example and numerical evidence 

In this section, we test the effectiveness of the MISC approximation on some instances of the general 
elliptic equation (3) in Example 1; more precisely, we consider a problem with one physical dimension (d = 1) 
as well as a more challenging problem with three dimensions (d = 3); in both cases, we set 23 = [0, l] d , 
fF(cc) = 1. As for the random diffusion coefficient, we set 

N 

a(x,y) = e lN( - x ' v \ Jn^x, y) = ^ K^n(x)y n , (28) 

n =1 

where y n are uniform random variables over T n = [—1,1], X n = v / 3exp(— n) and take f> n to be a tensorization 
of trigonometric functions. More precisely, we define the function 

sin irx^ if n is even 

cos ^——— itx^ if n is odd 

and set V’n(^) = <!>n(x) if d = 1. If d = 3, we take ipn( x ) = 4>i( n )( x i)4 , j(n)( x 2)<t , k(n){ x 3) f° r some indices 
i(n),j(n),k(n) detailed in Table 1. Observe that the boundedness of the supports of the random variables 
y n guarantees the existence of the two bounding constants in equation (4), a m j n and a max , that in turn 
assures the well posedness of the problem. Finally, the quantity of interest is defined as 

F (y) = J^ u {x,y)Q(x)dx, Q(x)= exp ^ 2cr 2 °^ 2 ) ( 29 ) 

with a = 0.16 and locations Xq = 0.3 for d = 1 and Xq = [0.3,0.2, 0.6] for d = 3. We also make the 
choice ho = 1/3 in Assumption 2. With these values and using the coarsest discretization, ho = 1/3, in 
all dimensions, the coefficient of variation of the quantity of interest can be approximated to be between 
90% and 100% depending on the number of dimensions, d, and the number of random variables, N, that we 
consider below. 

5.1. Verifying bounds on work and error contributions 

In this subsection we discuss the validity of Assumptions 2 to 4, upon which the MISC convergence 
theorem is based. To this end, we analyze separately the properties of the deterministic solver and of the 
collocation method applied to the problem just introduced. 



11 
















n 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

i{n) 

1 

2 

1 

1 

3 

2 

2 

1 

1 

1 

j(n) 

1 

1 

2 

1 

1 

2 

1 

3 

2 

1 

k(n ) 

1 

1 

1 

2 

1 

1 

2 

1 

2 

3 


Table 1: Included functions for d = 3 in (28). Here ip n {x) = 4>i( n ){ x i) ( t > j(n)( x 2)4 > k(n){ x 3)- 


Deterministic solver. The deterministic solver considered in this work consists of a tensorized finite dif¬ 
ference solver, with the grid size along each direction, x\,... ,Xd , defined by hi^ ai = h 0 2~ ai and no other 
discretization parameters are considered: therefore, D = d, Assumption 2 is satisfied, and, due to the 
Dirichlet boundary conditions prescribed for it, the overall number of degrees of freedom of the correspond¬ 
ing finite difference solution is nf=i(si--i)<nf=i( h 1 ^ . The associated linear system is solved with 
the GMRES method. We have numerically fitted the parameters, id and CgmreSj in the model: 

D 

Work[F a ] < Ggmres 

i=i 

for each individual tensor grid solve and found that d = 1 gives a good fit in our numerical experiments. 
From this we can recover the rates {7t}fLi and the constant C^ k in (19) with the following argument: 
since computing A dGt [F“] requires up to 2 D solver calls, each on a different grid (cf. equation (11)), we 
have 


Alb" = Work[A det [F a ]] = Work [F a _j] 

J'6{0,1} D 

D -i) 

< Cqmres ^2 IT 
16{0,1} D *= 1 

= Ggmres (^o2 _ai ) ^ ^2 H2 _J,,? 

\*=i / je{o,i} D i=i 

D 

= Cgmres(1 + ‘ 2 ~' } ) D H(/ii,Q i ) _,? , 

2=1 

i.e. , bound (19) is verified with 7 * = d,Mi = 1,..., D and C d ®* k = Cgmres(1 + 2~' a ) D . Hence, the sum of 
costs of the solver calls is proportional to the cost of the call on the finest grid. 

Concerning the error contribution AF det , we observe numerically that bound (20) holds true in practice 
with 77 = 2, * = 1,..., D, due to the fact that a £ C'°°(23) for p-almost every y £ T, T £ C°°(T>) and 
the function Q appearing in the quantity of interest (29) is also infinitely differentiable, confined in a small 
region inside the domain and zero up to machine precision on the boundary. In more detail, assuming for 
a moment that Assumption 3 is valid (we will numerically verify it later in this section), in Figure 2(a) we 
show the value of AE a p = AF det A£)g toc for fixed /3 = 1 and variable a = ja + l,j = 1, 2,..., as well as the 
corresponding value of the bound (20) for AF? det . The line obtained by choosing a. = [1, 0, 0] confirms that 
the size of A£? det indeed decreases exponentially fast with respect to au, and by fitting the computed values 
of AF det we obtain that the convergence rate is fj = 2, as previously mentioned; analogous conclusions can 
be obtained for 02 and 0:3 by setting a = [0, 1, 0] (shown in Figure 2(a)) and a. = [0, 0, 1] (not shown). 
Most importantly, confirmation of the product structure of A£' det can be obtained by observing, e.g., the 
decay of Afor ci: = [ 1 , 1 , 0 ] and a. = [ 1 , 1 , 1 ]. 

Stochastic discretization. The interpolation over the parameter space is based on the tensorized Lagrangian 
interpolation technique with Clenshaw-Curtis points explained in Section 2.1, cf. eqs. (5) and (6). In 
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(a) For fixed (3 = 1 and variable ex = kex + 1. 



(b) For fixed ex = 1 and variable (3 = k(3 + 1. 


Figure 2: Verifying the validity of the bound (23b) for the value of \AE a ^\ for the test case with D = 3 
and N = 5. The dashed lines are based on the model in (23b) with f j = 2 for all i = 1,2,3 and gj as in 
Table 2 for j = 1,..., 5. The solid lines are based on computed values. 


particular, due to the nestedness of the Clenshaw-Curtis points, adding the operator A stoc [F a!| g] to the 
MISC estimator will require AWp oc = IlyLi (m(/3j) — m([3j — 1)) new collocation points, which in view of 
equation (6) can be bounded as 

f 1 if Pj = 1 

m(/3j) - m(Pj - 1) = < 2 if = 2 < Vj = 1,2, ..., 

if/3,- >2, 

provided that the set X is downward closed: Assumption 2 and bound (21) in Assumption 4 are thus verified. 
Observe that the nestedness of the Clenshaw-Curtis knots is a key property here: indeed, if the nodes are 
not nested AVF| toc is not uniquely defined, i.e., it depends on the set X to which A stoc [F„ )/3 ] is being added, 
see, e.g., [34, Example 1 in Section 3]. 

Finally, to discuss the validity of bound (22) for A Ep oc , we rely on the theory developed in our previous 
works [35, 34]. We begin by introducing the Chebyshev polynomials of the first kind 'Fq(y) for y £ [—1,1], 
which are defined by the relation 

v k g (cos(0)) = cos (qO), 0 < 9 < 7r, q £ N. 

Then, for any multi-index q £ N n , we consider the A-variate Chebyshev polynomials tF q (y) = JI^Li (2 In) 

and introduce the spectral expansion of / : [—1,1]^ —> R. over {tFq}q e N w i 

r N 1 

f(y) = fg = / f(y)^ q (y) II a 2 

3r „ = i V 1 - Vl 

Next, given any £ 1 , £ 2 , • • •, £jv > 1 we introduce the Bernstein polyellipse £^ = 11^=1 , where 

denotes the ellipses in the complex plane defined as 

£ n,t„ = ^z n £ C: (z) < + ^ n cos3m (z) < ^ n sinift, (j) £ [0, 27t)|, 

and recall the following lemma (see [34, Lemma 2] for a proof). 
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Lemma 1. Let f : [—1, l] w —» R, and assume that there exist £ 1 ,^ 2 , • ■ •, £/v > 1 such that f admits a complex 
continuation f* : C N —> R holomorphic in the Bernstein polyellipse with sup zg£?i \f*(z)\ < B 

and B = ■ • ■ , £n) < 00. Then f admits a Chebyshev expansion that converges in C'°([— 1, l] w ), and 

whose coefficients f q are such that 


N 

\f q \ < Ccheb(q) e -3 " 9 ", g* n = log£ n 

n= 1 


(30) 


with Ccheb(g) = 2 M°B, where ||q||o denotes the number of non-zero elements of q. 

The following lemma then shows that the region of analyticity of F(y) indeed contains a Bernstein 
ellipse, so that a decay of exponential type can be expected for its Chebyshev coefficients. 

Lemma 2. The quantity of interest F{y ) = 0[u(-, y)\ is analytic in a Bernstein polyellipse with parameters 
f,n=T n + -\An + 1 , f or any r n < . 

Proof. Equation (3) can be extended in the complex domain by replacing y with z £ C N , and is analytic 
in the set £ = {z £ C N : $te[a(x, z)\ > 0 } , see, e.g., [ 6 ]. By writing z n = b n + ic n , we have 


i( X, z) = exp ^2 ZnKipnix)^ = exp ( y, b n X n ip n {x exp ic n X n ip n (x) 
= exp y b n X n ip n {x) cos PT c n X n ip n (x) + i sin ^ c n X n %f n (x) 

\ n / L \ n ) \ n / 


so that the region £ can be rewritten as 


£ = |z = 6 + ic £ C N : cos c n X n if n (x)^ > 0, Vx £ 231 . 


Such a region includes the smaller region 


£ 2 = { z = b + ic£ C N : 


y ' c n x n if r 


L°°m 


7r 

<2 


which in turn includes 


£3 — < z — b + ic £ C N : y> ' X n \c n \ < — 


where the last equality is due to the fact that, by construction, \\ifn\\ l°°(%) = 1> °f- equation (28). Next we 
let r n = 2 A C A and define the following subregion of £3: 

£4 = \z = b + ic£ C N : IcJ < T n \ C £3. 


£4 is actually a polystrip in the complex plain that it in turn contains the Bernstein ellipse with parameters 
t; n such that 

^ n ^ n = Tn => £,n - 1 - 2 T n £ n = 0 => = T n + + 1 

in which u{x , y) is analytic. Finally, the quantity of interest, F = 0[it], is also analytic in the same Bernstein 
polyellipse due to the linearity of the operator 0 . □ 
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9i 

92 

93 

94 

95 

96 

97 

5s 

59 

5io 

2.4855 

2.8174 

4.5044 

4.1938 

4.7459 

6.8444 

7.1513 

7.8622 

8.6584 

9.4545 


Table 2: Values of rates g for the test cases considered. 


Remark 9. Incidentally, we remark that the choice of r n considered in Lemma 2 degenerates for N —» 
oo. In this case, if we know that [C^o(^™IIVvlU°°(:B)) p < 00 f or some V < 1> then we could set r n = 
§(A„||^ ti ||l~(® ) ) p , which does not depend on N. 

Lemma 3. 

A E% oc < C E e~ 3lm(p n -l) M m{l3) ( 31 ) 

holds, where C E = 4, N B JX ^ =1 1 _ e 1 - 3 * • B as in Lemma 1, g* = log£„ with £ n as in Lemma 2, and 
has been defined in equation (7). 

Proof. Combining Lemmas 1 and 2, we obtain that the Chebyshev coefficients of F can be bounded as 

N 

\F q \<C Cheb (q)l[e-^, 

n—1 

with g* = log£ ra = log(r„ + \Jt% + 1) and r n as in Lemma 2. Then, the result can be obtained 
the same argument of [34, Lemma 5]. 

To conclude, we first observe that grows logarithmically with respect to m(/3), see eq. 

is asymptotically negligible in the estimate above, i.e. we can write 

N 

A£| oc < C E 2 (e) H 

n =1 

for an arbitrary e E > 0 and with C E2 (e E ) > C E , and furthermore that the definition of m{i) in ( 6 ) implies 
that m(i — 1) > . We can finally write 

N 

A£| oc < Ce 2 (e) JJ e -gn( 1 -«»)' m(/ 8 g ) ~ 1 

n =1 

with C|J:°o r = C(e) and g n = ^(1 — e E ). The latter bound actually shows that bound (22) in 

Assumption 4 is valid for the test we are considering. Finally, we point out that in practice we work with 
the expression (23b), whose rates g n are actually better estimated numerically, using the same procedure 
used to obtain the deterministic rates fj = 2 : we choose a sufficiently fine spatial resolution level a , consider 
a variable (3 = j(3 + 1 and fit the (simplified) model A Ep oc < C e -^ " 2,3 " . The values obtained are 
reported in Table 2, and they are found to be equal for the case d = 1 and d = 3 (see also [52, 35, 8 ]). To 
make sure that the estimated value of g n does not depend on the spatial discretization, one could repeat the 
procedure for a few different values of a and verify that the estimate is robust with respect to the spatial 
discretization: we note, however, that a rough estimate of g n will also be sufficient, since the convergence 
of the method is in practice dictated by the deterministic solver, as we have already discussed in Remark 7. 
Figure 2(b) then shows the validity of the bound A Ep oc < C n ^ =1 e -9 " 2 " comparing for fixed a = 1 and 
(3 = j(3 + \ the value of A£^[ et AEp oc and the corresponding estimate. 

Stochastic-deterministic product structure. We conclude this section by verifying Assumption 3, i.e., the fact 
that the error contribution can be factorized as A E a p = A£^ ot AE^ oc and that an analogous decomposition 
holds for AW a , p . While the latter is trivial, to verify the former we employ the same strategy used to verify 
the models for Aand A Ep oc , this time letting both a and (3 change for every point, i.e., a = jot + 1 
and (3 = j/3o + 1. Figure 3 shows the comparison between the computed value of A E a ,p and their estimated 
counterpart and confirms the validity of the product structure assumption. 


N 

^stoc TT -g„m(t3„) 
error 1 5 

n—1 


following 

□ 

( 8 ), so it 
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Figure 3: Comparison of \AE a ^p\ for /3 = j/3 + 1 and a = jet + 1 for the test case with D = 3 and N = 5. 
The dashed lines are based on the model in (23b) with r i = 2 for all i = 1,2,3 and gj as in Table 2 for 
j = 1,..., 5. The solid lines are based on computed values. 

5.2. Test setup 

In our numerical tests, we compare MISC with the methods listed below. For each of them we show (for 
both test cases considered) plots of the convergence of the error in the computation of E[.F] with respect to 
the computational work, taking as a reference value the result obtained using a well-resolved MISC solution. 
To avoid discrepancies in running time due to implementation details, the computational work is estimated 
in terms of the total number of degrees of freedom, i.e., using (14) and (23a). The names used here for the 
methods are also used in the legends of the figures showing the convergence plots. 

“a-priori” MISC refers to the MISC method with index set I defined by (17), where AW a ^p and AE a p 
are taken to equal their upper bounds in (23a) and (23b), respectively. The resulting set is explicitly 
written in (24). The convergence rate of this set is predicted by Theorem 1, cf. Remark 6. Note that 
we do not need to determine the value of the constants Cwork and C er ror since they can be absorbed 
in the parameter e in (17). 

“a-posteriori” MISC refers to the MISC method with index set 1 defined by (17), where AW a ,p is taken 
to equal its upper bound in (23a), and A E a p is instead computed explicitly as |A[F Qi/ 3 ]|. Notice that 
this method is not practical since the cost of constructing set T would dominate the cost of the MISC 
estimator by far. However, this method would produce the best possible convergence and serve as a 
benchmark for both “a-priori” MISC and the bound (23b). 

MLSC (only in the case d > 1) refers to the Multilevel Stochastic Collocation obtained by setting a\ = 
... = Qfl (i.e. considering the mesh-size as the only discretization parameter), as already mentioned 
in Remark 5; we recall this is not exactly the MLSC method that was implemented in [24], see again 
Remark 5. Just as with MISC, we consider both the “a-priori” and “a-posteriori” version of MLSC, 
where A E a p is taken to be equal to its upper (23b) in the former case and assessed by direct numerical 
evaluation in the latter case. 

SCC refers to the “Sparse Composite Collocation method” in Remark 5, see equation (18). 

MIMC refers to the Multi-Index Monte Carlo method as detailed in [29], for which the complexity 
O can be estimated for the test case at hand and as long as d < 4. 

SGSC refers to the quasi-optimal Sparse Grids Stochastic Collocation (SGSC) with fixed spatial discretiza¬ 
tion as proposed in [35, 34]. To determine the needed spatial discretization for a given work and for a 
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fair comparison against MISC, we actually compute the convergence curves of SGSC for all relevant 
levels of spatial discretizations and then show in the plots only the lower envelope of the corresponding 
convergence curves, ignoring the possible spurious reductions of error that might happen due to non- 
asymptotic, unpredictable cancellations, cf. Figure 4. In this way, we ensure that the error shown for 
such “single-level methods” has been obtained with the smallest computational error possible. Again, 
this is not a computationally practical method but is taken as a reference for what a sparse grids 
Stochastic Collocation method with optimal balancing of the space and stochastic discretization errors 
could achieve. 



Figure 4: Envelope of SGSC convergence curves for the test case with d = 3 and N = 10. 


5.3. Implementation details 

To implement MISC, we need two components: 

1. Given a profit level parameter, e, we build the quasi-optimal set 1 based on (17), (23a) and (23b). 
One method to achieve this is to exploit the fact that this set is downward closed and use the following 
recursive algorithm. 

FUNCTION BuildSet(epsilon, multiindex) 

FOR i = 1 to (D+N) 

IF Profit(multiindex + e_i) > epsilon 
THEN 

ADD multilndex+e_i to FinalSet 

CALL BuildSet(epsilon, multilndex+e_i) 

END IF 
END FOR 
END FUNCTION 


2. Given the set, I(L), we evaluate (13). Here we have two choices: 

• Evaluate the individual terms A[F a p] for every at, (3 £l. To do so, we use the operator defined 
in (9) along each stochastic and spatial direction. By storing the values of these terms, we can 
evaluate the MISC with different index sets (contained in I{L)), which might be required to 
test the convergence of the MISC method. Moreover, this implementation is suitable for adaptive 
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methods that expand the index set based on some criteria and reevaluate the MISC estimator. On 
the other hand, this implementation has a computational overhead since most computed values of 
Fa,/3 will actually not contribute to the final value of the estimator. However, this computational 
overhead is only a fraction of the minimum time required to evaluate the estimator. 

• Use the combination form of (13) and only compute the terms that have c a p 7 ^ 0. This would 
remove the overhead of computing terms that make zero contribution to the estimator. This 
implementation is more efficient but less flexible as we cannot evaluate the estimator on sets 
contained in I(L) or build the set adaptively. 

5-4■ Test with D = 1 

Here we consider three different numbers of stochastic variables, namely N = 1, 5,10. Results are shown 
in Figure 5. As expected, a-posteriori MISC shows the best convergence, with a-priori MISC being slightly 
worse and the single level methods following. Finally, we verify the accuracy of the estimated asymptotic 
convergence rate provided by Theorem 1: in this case, C = r J = D ^ 2 = 2 and 3 = 1 holds. Hence, the 
predicted convergence rate is W _< >(log fU)^ + b( 3_1 ) = W~ 2 , which appears to be in good agreement with 
the experimental convergence rate. 

5.5. Test with D = 3 

In this case, we obtain the convergence curves shown in Figure 6 , where the Multilevel Stochastic 
Collocation method has also been included. The hierarchy between the methods is in agreement with the 
case d = 1, with the Multilevel Stochastic Collocation being comparable or slightly better than single level 
methods, but worse than the MISC approaches as expected. 

Concerning the accuracy of the theoretical estimate: since for this test r 1 = r 2 = r 2 = 2 and 
7 i = 72 = 73 = 1, C = 2 still holds, while this time 3 = 3; hence, the predicted convergence rate is 
IU -< ’(logfU)^ +1 ^ 3-1 ) = W _ 2 (logW) 6 . The plots suggest that the theoretical estimates might be slightly 
too optimistic when N increases but it is important to recall that Theorem 1 gives only an asymptotic 
result, and the plot could be negatively influenced by pre-asymptotic effects. Observe also that in this case 
there are a few data points where a-posteriori MISC is not better than a-priori MISC; this observation can 
be ascribed to the fact that a-posteriori MISC is optimal only with respect to the upper bound in (15). In 
other words, a-posteriori MISC selects the contributions according to the absolute value of the contributions 
but then the MISC estimator is computed by summing signed contributions. Hence, cancellations between 
contributions with similar sizes and opposite signs will occur. 

Finally, we remark that, in our calculations, MLSC and SGSC were not able to achieve very small errors, 
unlike MISC. This is due to a limitation in the linear solver we are using that allows systems with only up 
to 2 17 degrees of freedom (around 1GB of memory) to be solved. These “single-level” methods hit that limit 
sooner than MISC since they entail solving a very large system that comes from isotropically discretizing 
all three spatial dimensions. 

6. Conclusions 

In this work, we have proposed MISC, a combination technique method to solve UQ problems, optimizing 
both the deterministic and stochastic resolution levels simultaneously to minimize the computational cost. 
A distinctive feature of MISC is that its construction is based on the notion of profit of the mixed differences 
composing it, rather than balancing the total error contributions arising from the deterministic and stochastic 
components. We have detailed a complexity analysis and derived a convergence theorem showing that in 
certain cases the convergence of the method is essentially dictated by the convergence properties of the 
deterministic solver. We have then verified the effectiveness of the method proposed on a couple of numerical 
test cases, comparing its performance with other methods available in the literature. The results obtained 
are encouraging, as they suggest that the proposed methodology is more effective than the other methods 
considered here. The theoretical results have been also found to be consistent with the numerical results to 
a satisfactory extent. 
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Figure 5: Results for test D = 1, case N = 1 (top), N = 5 (center) and N = 10 (bottom). 
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Figure 6: Results for test D = 3, case N = 1 (top), N = 5 (center) and N = 10 (bottom). 
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As a final remark, we observe that the methodology presented here is not limited to the spatial or temporal 
discretization parameters of the deterministic problem, but could also be applied to other discretization 
parameters, such as smoothing parameters or artificial viscosities. 
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Appendix A. Proof of Theorem 1 

The following technical lemmas are needed in the convergence proof. 

Lemma 4. For x G (1, oo) D , define \x\ = (\xi\ )^ =1 . For any f : (1, oo) D -4 R and g : (1, oo) D —> R+, 


Y 9 (a) = 


(«eN° : /(«)<o} J{*e(i,oo)° : /(L*J)<o} 


g([x\) dx 


holds. Moreover, if g and f are increasing, then 

Y 9(a) < 


{ace : /(tt)<0} 


{x^(l,oo) D : f(x— 1)<0} 


g(x) dx, 


and if g and f are decreasing, then 

Y 9(a) < 


{ckGNJ : /(c*)< 0 } 


(ccG^oo) 0 : f(x)< 0} 


g(x — 1) dx. 


Proof. We have 


Y 9(a) = Y 9(a) 

{a£N® : /(a)< 0 > {a£N“ : /(a)< 0 } 


da; 


Y / 3 (|a + ®J)da; 

, JxG[ 0 , 1 ]° 


{a£N® : /(a)< 0 > 


J {x£(1,oo) d : /(La=J)<0} 

Combining these inequalities with x — 1 < < x finishes the proof. 

Lemma 5. Assume a G R£, b G R£ and L > |a|. Then, 

/ d \ / D 


g(\x\) dx. 


Y 

{a;£N+ : YiLi o,ie b i x i-\-biXi>L} 


exp (- y aiebixi ) - (n cxp ^ a ^ ) exp (~ l ) ( l +i) 2D+i . 


Proof. Define the set 


V = 




□ 
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and define [y\ = {\_Vi\)f =v Then 


D 


E ex p -E 


n biXi 


{ajGNj : EiLi ^ie biX -\-biXi>L} 


ai e | = Y exp (—a • y) 

i=l / {yeP : a-y+l log(y)|>£} 

< Y exp(-a|yj) 

{yeP : a-(Lj/J+l) + | log(Ll/J+l)|>L} 

< E 

+ ’ ~ 1 


< 


J {ye(l,oo) D : a-y+l log(y+l)|>L-|a|} 

Letting z t = a^i + log(y.j + 1) and p(zi) = yt < + then 

/ d \ 


exp {-a ■ y ) 

exp(—a • (y - 1)) dy. 


E 


exp 


-E 




biXi 


{aj£N+ : J2?=i a,ie b i x +biXi>L} 


i=1 


= exp(|a|) 
= exp(|a|) 


{j/e(l,c»)- D : a-y+| log(y+l)|>L-|a|} 


exp(—a -y- | log(y + 1)| + | log(y + 1)|) dy 


D 


' {«e0|^ 1 (ai+log(2),oo) : |«|>L—|a|} 


e xp+|z|)n( P ^' )+ 1 1 ) dz 

i=i \ a * + R5IHT/ 


D 


< exp(|o|) 


' {«e0^ 1 (ai+log(2),oo) : |z|>L— |a|} 


exp(-|«|)p ( j dz 


< 


/ D 

n 


exp(oi) 


_1 a i J J {ze®P =1 (cii+\og(2),oo) : \z\>L—\a\} 
D 


exp(—|z| + | log(z + a)|) dz 


n 


\*=1 

*(n 


exp (2 a,;) 
exp (2 a,;) 


/ exp(—1*| + |log(x)|) dec 

'(a:S®^ 1 (log(2),oo) : |x|>L} 

/ exp(—|z| + | log(z)|) dz. 

' {«G0^ 1 (O,oo) : |«|>L} 


Now let us prove, by induction on D, that we have 

^ exp(—|z| + | log(z)|) dz < exp(-L)(L + 1) 2£> ~ 1 


'{z£l® : |jz| >L} 


For D = 1, the inequality is a trivial equality that can be obtained with integration by parts. Assume the 
inequality is true for D and let us prove it for D + 1: 


l{ze R+ +1 : |z|>£} 


exp+|z| + log(z)) dz = / yexp(— y) / exp(—|er| + log(ec)) dee 

Jl J{xe R+} 


dy 


+ / yexp(-y) / exp(-|ec| + log(x)) da? dy 

J 0 J{xdzRP : \x\~>L—y} 


< exp {-L)(L + 1) 


+ [ V exp(—y) exp(—L + y)(L — y + 1) 2D 1 dy 
Jo 
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< exp(-L) (L + 1 + L 2 {L + l) 2D ~ l ) 

< exp(-L) (L + 1) (1 + L{L + l)^" 1 ) 

< exp(—L)(L + l ) 2 ( D+1 ) -1 . 


Finally, substituting back, we get the result. 


□ 


Definition 1. Given a £ and A > 0, let n(a, A) denote the number of occurrences of A in a, 

n(a, A) = #{i = 1,..., d : = A}. 

Lemma 6. Assume k £ N, a £ R+,L > |a|. Then, the following bounds hold true: 

[ exp (—a • x) dx< 55 D {a ) exp(- rnin(a)L)L n ( a ’ min ( a »- 1 

where 55 d{o) is a positive constant independent of L. 

Proof. See [29, Lemma B.3] for a proof of the inequality and the value of 55 d(ci). Moreover, a proof of a 
consistent equality for the case a = 1 can be found in [51, Proposition 2.3]. □ 

Lemma 7. Assume k eN, a £ R+,L > |a|. Then, the following bound holds: 

exp (a ■ x) (L — |cc|) fe dx < 2 l£>( 0 ', k) exp(max(a)L)L I '^ a,max( '“' ) ' 1 ” 1 , 


where 


' : \x\<L} 


21 D (a,k) = 


k\ 


(n(a,max(a)) — l)!max(a ) fe+1 


n 


i 


iai<max(a) 


rnax(a) — at 


Proof. Without loss of generality, assume that a,; > a i+ 1 for all i = 1... D, such that oq = max(a). We 
prove the result by induction on D. For D = 1, we have 

f exp {ax) {L - x) k dx = (exp {aL) - j 


i =0 


< 


k\ exp {aL) 


1 


Next, assume that the result is valid for a given D > 1 and a £ where aq > Gq+i for alii = 1... D, 
such that a\ = max(a). Let b < a\ and define a new vector a = {a, b) £ R^ +1 . We have 


'{(*,!/)eR ° +1 : y+\x\<L} 


exp {by + a - x) {L — y — |a|) dy da; 


= / exp {by) / exp (a • x) {L — y — |a;|) fc da; dy 

Jo J {a : \x\<L—y} 


< 2l D (a,fc)exp(eqL) [ exp {{b - ai)y){L - y) n[a ' ai) 1 dy. 
Jo 

We distinguish between two cases: 
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1 . b < a\ then n(a, ai) = n(a, oi) and 


pL nOO 

/ exp (-( ai - b)y)(L - y)"M-l d y < L n ^~ l / exp (-( 0l - b)y) d y 
Jo Jo 


< £ n ( a >°i) _1 


ai — b ' 


and in this case 


21 D(a,k) 


k\ 


a\ — b (n(a, ai) — l)!aj +1 
= 21 d+ i (a,k). 


n 

\di<ai 


a\ — ai ) \ai - b 


2 . b = a± then n(a, ai) = n(a, ai) + 1 and 

pL 


(L-y) 


n(a,ai)-l ^ _ 


j^n(a,ai) j^n{a,a\)— 1 

n(a, m) n(a, ai) — 1 ’ 


and again 


n (a, ai) — 1 


21 D {a,k) = 


k\ 


n(a, ai) — 1 (n(a, a\) — l)!a ^ +1 

= 2 l_D+i(a, k). 


n 

\ai<a\ 


a\ — ai 


□ 


Theorem 1 (MISC computational complexity). Under Assumptions 2 to f, the bounds for the factors 
appearing in Assumption 3 can be equivalently rewritten as 


AW r a,fl<Cw 0 rk e E « 7 ,B< e , l' J l ) 

AE at p < C erioI e~ ^.= 1 e- ^ 9 > exp(5ft ) , 
with 7 i = t i log 2, r,; = r* log 2, S = log 2 and gi = 'giC rn y ow . Define the following set 

T*(L) = I [a,/3] € n° +N : ^(r 4 + 7 ;H + ^(<% +^e 5ft ) < i] with L e R+. 

I 2=1 2=1 J 

Then there exists a constant Gw such that, for any W max satisfying 

w max > Civ exp (x), 

and choosing L as 


L = L(W max ) = - ( log 
X 


w E 


• w 


- (3 - 1) log - log 
X 




'W 


(23a) 

(23b) 


(24) 


(25) 


(26) 


with 


-(= 


7i 


7l+ri , • • • , 7o 7 + r ^ ) > X = max(E) , C = min^i,...^ ^ and i = #{i = 1,... D : ^ = (}, MISC 
estimator Mx* (£(w max )) satisfies 


Work[Mx*(i(W max ))] < Umax, 

Error jj] 

hmsup -:---TTTTu-TT 

wwr°° w m L (log (W max )) (C+ )(3 


= 6 e < oo. 


(27a) 

(27b) 
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Proof. The bounds (23a) and (23b) can be obtained by elementary algebraic operations combining Assump¬ 
tions 2 and 4; for instance, 


D 


D 


D 


D 


AWT < c£ k nv* = c wor k n( /i ° 2_ “ , )^' : = II 2 ^ = Ip ia,los2 


iV 


i=l 

N 


a w|° c < n 2/3,1 = c wo° c k n 


,Pn log 2 


n=l 


n—1 


from which (23a) follows by setting Cwork = C^®* k /i 0 ^^work- The P r °of is then divided into two steps. 

Step 1: Work Estimate. Observe that H* = ^ r . < 1 for all i = 1, ... D and that 3 = n(E, y). Thanks to 
equations (14) and (23a), and using Lemma 4, the total work satisfies 


Work [X* (L)] = J2 AW <*,P 

(a,/9)er*(£) 

A O\vork ^ ^ 

{(a,/3)eN° +JV : Ef = i(n+7b«,+Ef=i S^+gje^JKL} 

— f\vork 


/{(«,/3)e(i,oo)°+« : E,=i0i+7i)O,-i)+Ef = 1 

Next, let /3j = gje s ^ j ~^ and al = {ji + ji)(ai — 1). We have 


ex p [fT/liUi + <5|/3|^ 

exp y] 7 iQfi + 5\(3\ do; d/3. 


Work [I* (I/)] <C work 




exp(7i) 

n + 7 i 


/{(«,/3)eR?x( 


Dropping the over-line notation and defining L 


exp (E • a) dec d/3. 

: I“l + I0| + | l°g/3|<i+l logsl} 

= X + | log fir | and Cw,i to be the constant factor, we obtain 


Work[X*(X)] < e w ,i 


= &w,i [ 

h 


'{(a,/3)6R^x(®Jk 1 (gj,oo)) : |t»| + |/3|+| log/3|<L} 


exp (E ■ a) da d/3 


{/3e®|L 1 (sj, 00 ) : 1/31+1 log/3|<i} J{a£ R£ : \a\<L~\/3\-\ log/3|} 


exp (E ■ a) da d/3 


< Cw.iSId (2,0) / 

3{ Pe®f =1 (gj,oo) : |/3| + | log/3 |<l} 

Define Cw ,2 = Cw,i21d (3, 0) exp(x’X), then 


exp (x (l - \(3\ - | log/31 ) ) (x - |/3| - | log/3|) d/3 


Work[X*(X)] < C w ,2 
Gw,2 ( 


< 


{/3e®jX 1 (fl ; ,',cx>) : |/3| + | log/3|<L} 

L-\g \-I logfipl) 3 


exp/— 


(—X (|/3| + | log /31)) (x — |/3| — | log /31) 3 1 d/3 


'{/3e®y =1 (gj,co) : |3[+| log/3|<i} 

Since x > 0, the previous integral is bounded for all L and we have 

Work[X*(I/)] < C w exp(xL) (X - |g|) 3_1 < C w exp(x'X)X 3_1 , 


exp(-x(l/3| + |log/3|)) d/3. 
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where 


N 


C\V — Cwork I ] J[ 




g : > k) g 2 

Substituting (26) yields 


/ D 

n 


exp fa) 
ft +7* 


21 d (3,0) 


' {/3e®fa(sj ,oo)} 


exp(-x(|/3| + |log/3|)) d/3. 


Work [I*(L)] < kfa 


(d - 1 ) log 


1 - 




3-1 


V 


lo S (^r) 




From here it is easy to see that if (25) is satisfied, then (27a) follows. 

Step 2: Error Estimate. Thanks to equations (15) and (23b), the total error satisfies 
Error [I* (L)} < E A E a ^ 


D 


N 


— Herror ^ ^ 

{(a,/3)efa +JV : E”i(r.+7.)a.+E"i%+»e^>L} 

D N 


exp _ E riai - E g i e 




»=1 l=i 


f Vrru ;■ ^ ^ 

{(a,/3)eN“ +JV : £? = i(r-,+ 7 ,K>£} 


ex P -E^< 


,«& 


i=! 1=1 


D 


AT 


c error e E 

{aeN° : £f = i(n+7«)«*<£} {/36N« : £f =1 ift+Sje^^i-Ef^tri+Tila,} 


exp - E ri0ii ~ E g i e 


Mi 


i =i l=i 


Looking at the first term, let 77 .j = fa r < 1 and r/ = fa ) i=1 and note that 3 = # {/ = 1... D : rji = mmfa)}. 
Then 


E 

{(a,/3)eN^ +JV : Ef=,(ri+ 7 <)«i>i} 


D 


N 


exp -E r ^ - E** 


,«&• 


i=i 7=1 


A/ - 


/ 


= E exp -E^ e 




V/3GFft 


1=1 


E 


y{aeN° : E(Li( r i+7i)«i>i} 


exp - E ri0li 

V i=l 1 


£ e&1 l 


D 


{aG(l,oo) D : EiLi ( r i +li)<Xi > L } 


exp -E^-i) da 


D 


expfa 


D 


= Ce,i { JJ . . 

\j = l r i ' li J J {cce®(7 1 (ri+7i,oo) : \x\>L} 

< e E , 2 exp (- mm(r])L) L 3_1 , 


exp - E 




dx 


where 


D 


fa, 2 = 23d(t 7) j n 


expfa) 
ft + 7, 


AT 


E exp -E& e4ft 


/3GN" 


1=1 
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For the second term, letting H = L — YliLi( r i + 7 we can bound the sum using Lemma 5: 


E 


N 


N 


exp ( ] < ( II 

3= 1 


exp (2 gj) 


exp (—H)(H + 1) 


2N— 1 


{/ 3 GN« : Ef =1 Sflj+gje^^H} 

Defining C_e ,3 = IljLi ex P(^9j)9j 2 and substituting back 


= 1 


D 


E exp ~ E Ti<Xi 


E 


N 


exp E^< 




{aeN° :Ef =1 (ri+70«»<i} V i=1 7 {/3£N« :Ef=i'5ft+We 5 ^>L-E"i(n+7i)«i} V j=1 


D 


D 


2N-1 


<Ge , 3 ^2 exp I -L + 7 jag j I L + 1 - Y/ r » + 

{q£H“ : ES=i(n+7i)«i<£} ' i=1 7 V *=1 > 

„ / D \ / D 


2JV-1 


= e 


_E,3 


{o;e(l,oo)r> : EiLi( r 'i+7i)L“iJ< Z '} 


exp -L + ^2 7, |_£bj L + 1 - + 7,) (a*) 


da 


D 


i= 1 

D 


2N-1 


< Qe,: 


{«G(l,oo)-° : EiLi( r i+7i)(Q!i-l)<^} 

f 

{aeRj : |x|<i} 


exp —1/ + ^7*07 L + 1 - ^(77 + 7»)(ai - 1) 


da 


< e E ,4exp((x- 1)L)L 3_1 , 


exp (S • £c) (L + 1 — |®|) 2Ar 1 dec 


where 


N 


D 


e E , 4 - (n 

<3 =1 


exp ( 2 5j) 1 / TT exp (7i) 




n 


\2 — 1 


Finally, noting that 


7 i + U 

X - 1 = -min(rj), 


sa D (s, 2 JV-i). 


we have the error estimate 

Error [X*(L)] < C e „ OT (C Ei2 + e Ej4 ) exp(-mm(r/)L)L 3_1 . 
Then, substituting L from (26) and evaluating the limit gives (27b). 


□ 


References 

[1] R. G. Ghanem, P. D. Spanos, Stochastic finite elements: a spectral approach, Springer-Verlag, New York, 1991. 

[2] O. P. Le Maitre, O. M. Knio, Spectral methods for uncertainty quantification, Scientific Computation, Springer, New 
York, 2010, with applications to computational fluid dynamics, doi: 10.1007/978-90-481-3520-2. 

[3] H. G. Matthies, A. Keese, Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations, 
Computer Methods in Applied Mechanics and Engineering 194 (12-16) (2005) 1295-1331. 

[4] R. A. Todor, C. Schwab, Convergence rates for sparse chaos approximations of elliptic problems with stochastic coefficients, 
IMA J Numer Anal 27 (2) (2007) 232-261. doi: 10.1093/imanum/drl025. 

[5] D. Xiu, G. Karniadakis, The Wiener-Askey polynomial chaos for stochastic differential equations, SIAM Journal on 
Scientific Computing 24 (2) (2002) 619-644. 

[6] I. Babuska, F. Nobile, R. Tempone, A stochastic collocation method for elliptic partial differential equations with random 
input data, SIAM Review 52 (2) (2010) 317-355. 

[7] H. J. Bungartz, M. Griebel, Sparse grids, Acta Numerica 13 (2004) 147-269. 

[8] F. Nobile, R. Tempone, C. Webster, An anisotropic sparse grid stochastic collocation method for partial differential 
equations with random input data, SIAM Journal on Numerical Analysis 46 (5) (2008) 2411-2442. 


27 






[9] D. Xiu, J. Hesthaven, High-order collocation methods for differential equations with random inputs, SIAM Journal on 
Scientific Computing 27 (3) (2005) 1118-1139. 

[10] B. Khoromskij, C. Schwab, Tensor-structured galerkin approximation of parametric and stochastic elliptic pdes, SIAM 
Journal on Scientific Computing 33 (1) (2011) 364-385. 

[11] B. Khoromskij, I. Oseledets, Quantics-tt collocation approximation of parameter-dependent and stochastic elliptic pdes, 
Computational Methods in Applied Mathematics 10 (4) (2010) 376-394. 

[12] A. Nouy, Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace 
problem and dedicated algorithms, Computer Methods in Applied Mechanics and Engineering 197 (51) (2008) 4718-4736. 

[13] J. Ballani, L. Grasedyck, Hierarchical Tensor Approximation of Output Quantities of Parameter-Dependent PDEs, 
SIAM/ASA Journal on Uncertainty Quantification 3 (1) (2015) 852-872. doi: 10.1137/140960980. 

[14] S. Boyaval, C. Le Bris, T. Lelievre, Y. Maday, N. Nguyen, A. Patera, Reduced basis techniques for stochastic problems, 
Archives of Computational Methods in Engineering 17 (4) (2010) 435-454. 

[15] P. Chen, A. Quarteroni, G. Rozza, Comparison between reduced basis and stochastic collocation methods for elliptic 
problems, Journal of Scientific Computing 59 (1) (2014) 187-216. doi: 10.1007/sl0915-013-9764-2. 

[16] A. Cohen, R. Devore, C. Schwab, Analytic regularity and polynomial approximation of parametric and stochastic elliptic 
PDE’S, Analysis and Applications 9 (1) (2011) 11-47. 

[17] S. Heinrich, Multilevel Monte Carlo methods, in: Large-Scale Scientific Computing, Vol. 2179 of Lecture Notes in Computer 
Science, Springer Berlin Heidelberg, 2001, pp. 58-67. 

[18] M. B. Giles, Multilevel Monte Carlo path simulation, Operations Research 56 (3) (2008) 607-617. 

[19] A. Barth, C. Schwab, N. Zollinger, Multi-level Monte Carlo finite element method for elliptic PDEs with stochastic 
coefficients, Numerische Mathematik 119 (1) (2011) 123-161. 

[20] A. Barth, A. Lang, C. Schwab, Multilevel Monte Carlo method for parabolic stochastic partial differential equations, BIT 
Numerical Mathematics 53 (1) (2013) 3—27. 

[21] J. Charrier, R. Scheichl, A. Teckentrup, Finite element error analysis of elliptic pdes with random coefficients and its 
application to multilevel Monte Carlo methods, SIAM Journal on Numerical Analysis 51 (1) (2013) 322-352. 

[22] K. Cliffe, M. Giles, R. Scheichl, A. Teckentrup, Multilevel Monte Carlo methods and applications to elliptic pdes with 
random coefficients, Computing and Visualization in Science 14 (1) (2011) 3-15. doi: 10.1007/s00791-011-0160-x. 

[23] S. Mishra, C. Schwab, J. Sukys, Multi-level Monte Carlo finite volume methods for nonlinear systems of conservation laws 
in multi-dimensions, Journal of Computational Physics 231 (8) (2012) 3365—3388. 

[24] A. Teckentrup, P. Jantsch, C. G. Webster, M. Gunzburger, A Multilevel Stochastic Collocation Method for Partial 
Differential Equations with Random Input Data, SIAM/AS A Journal on Uncertainty Quantification 3 (1) (2015) 1046- 
1074. doi:10.1137/140969002. 

[25] H. W. van Wyk, Multilevel sparse grid methods for elliptic partial differential equations with random coefficients, arXiv 
preprint arXiv: 1404.0963 (2014). 

[26] H. Harbrecht, M. Peters, M. Siebenmorgen, On multilevel quadrature for elliptic stochastic partial differential equations, 
in: Sparse Grids and Applications, Vol. 88 of Lecture Notes in Computational Science and Engineering, Springer, 2013, 
pp. 161-179. 

[27] F. Y. Kuo, C. Schwab, I. Sloan, Multi-level Quasi-Monte Carlo Finite Element Methods for a Class of Elliptic PDEs with 
Random Coefficients, Foundations of Computational Mathematics 15 (2) (2015) 411-449. 

[28] F. Nobile, F. Tesei, A Multi Level Monte Carlo method with control variate for elliptic PDEs with log-normal co¬ 
efficients, Stochastic Partial Differential Equations: Analysis and Computations 3 (3) (2015) 398-444. doi: 10.1007/ 
S40072-015-0055-9. 

[29] A.-L. Haji-Ali, F. Nobile, R. Tempone, Multi-index Monte Carlo: when sparsity meets sampling, Numerische Mathematik 
132 (2015) 767-806. doi:10.1007/s00211-015-0734-5. 

[30] H. J. Bungartz, M. Griebel, D. Roschke, C. Zenger, Pointwise convergence of the combination technique for the Laplace 
equation, East-West Journal of Numerical Mathematics 2 (1994) 21-45. 

[31] M. Griebel, M. Schneider, C. Zenger, A combination technique for the solution of sparse grid problems, in: P. de Groen, 
R. Beauwens (Eds.), Iterative Methods in Linear Algebra, IMACS, Elsevier, North Holland, 1992, pp. 263-281. 

[32] M. Hegland, J. Garcke, V. Challis, The combination technique and some generalisations, Linear Algebra and its Applica¬ 
tions 420 (2-3) (2007) 249-275. doi:10.1016/j.laa.2006.07.014. 

[33] M. Griebel, H. Harbrecht, On the convergence of the combination technique, in: J. Garcke, D. Pfliiger (Eds.), Sparse 
Grids and Applications - Munich 2012, Vol. 97 of Lecture Notes in Computational Science and Engineering, Springer 
International Publishing, 2014, pp. 55-74. doi: 10.1007/978-3-319-04537-5_3. 

[34] F. Nobile, L. Tamellini, R. Tempone, Convergence of quasi-optimal sparse-grid approximation of Hilbert-space-valued 
functions: application to random elliptic PDEs, Numerische Mathematikdoi : 10.1007/s00211-015-0773-y. 

[35] J. Beck, F. Nobile, L. Tamellini, R. Tempone, On the optimal polynomial approximation of stochastic PDEs by Galerkin 
and collocation methods, Mathematical Models and Methods in Applied Sciences 22 (09) (2012) 1250023. 

[36] J. Beck, F. Nobile, L. Tamellini, R. Tempone, A Quasi-optimal Sparse Grids Procedure for Groundwater Flows, in: 
Spectral and High Order Methods for Partial Differential Equations - ICOSAHOM 2012, Vol. 95 of Lecture Notes in 
Computational Science and Engineering, Springer, 2014, pp. 1-16. 

[37] M. Griebel, S. Knapek, Optimized general sparse grid approximation spaces for operator equations, Mathematics of 
Computation 78 (268) (2009) 2223-2257. doi: 10.1090/S0025-5718-09-02248-0. 

[38] M. Bieri, A sparse composite collocation finite element method for elliptic SPDEs., SIAM Journal on Numerical Analysis 
49 (6) (2011) 2277-2301. 

[39] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, nurbs, exact geometry and mesh 


28 


refinement, Computer Methods in Applied Mechanics and Engineering 194 (39-41) (2005) 4135-4195. doi:10.1016/j. 
cma.2004.10.008. 

[40] W. J. Gordon, C. A. Hall, Construction of curvilinear co-ordinate systems and applications to mesh generation, Interna¬ 
tional Journal for Numerical Methods in Engineering 7 (4) (1973) 461-477. doi: 10.1002/nme. 1620070405. 

[41] A. Quarteroni, A. Valli, Domain Decomposition Methods for Partial Differential Equations, Numerical mathematics and 
scientific computation, Clarendon Press, 1999. 

[42] L. Trefethen, Approximation Theory and Approximation Practice, Society for Industrial and Applied Mathematics, 2013. 

[43] L. N. Trefethen, Is Gauss quadrature better than Clenshaw-Curtis?, SIAM Review 50 (1) (2008) 67-87. 

[44] A. Chkifa, On the lebesgue constant of leja sequences for the complex unit disk and of their real projection, Journal of 
Approximation Theory 166 (0) (2013) 176-200. 

[45] F. Nobile, L. Tamellini, R. Tempone, Comparison of Clenshaw-Curtis and Leja Quasi-Optimal Sparse Grids for the 
Approximation of Random PDEs, in: R. M. Kirby, M. Berzins, J. S. Hesthaven (Eds.), Spectral and High Order Methods 
for Partial Differential Equations - ICOSAHOM ’14, Vol. 106 of Lecture Notes in Computational Science and Engineering, 
Springer International Publishing, 2015, pp. 475-482. doi : 10.1007/978-3-319-19800-2_44. 

[46] A. Narayan, J. D. Jakeman, Adaptive Leja Sparse Grid Constructions for Stochastic Collocation and High-Dimensional 
Approximation, SIAM Journal on Scientific Computing 36 (6) (2014) A2952-A2983. 

[47] A. Genz, B. D. Keister, Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian 
weight, Journal of Computational and Applied Mathematics 71 (2) (1996) 299-309. 

[48] G. W. Wasilkowski, H. Wozniakowski, Explicit cost bounds of algorithms for multivariate tensor product problems, Journal 
of Complexity 11 (1) (1995) 1-56. 

[49] S. Martello, P. Toth, Knapsack problems: algorithms and computer implementations, Wiley-Interscience series in discrete 
mathematics and optimization, J. Wiley & Sons, 1990. 

[50] D. Dung, M. Griebel, Hyperbolic cross approximation in infinite dimensions, Journal of Complexitydoi: 10.1016/j .jco. 
2015.09.006. 

[51] M. Griebel, J. Oettershagen, On tensor product approximation of analytic functions, Journal of Approximation Theory 
(2016) In printdoi :10.1016/j.j at.2016.02.006. 

[52] J. Back, F. Nobile, L. Tamellini, R. Tempone, Stochastic spectral Galerkin and collocation methods for PDEs with random 
coefficients: a numerical comparison, in: Spectral and High Order Methods for Partial Differential Equations, Vol. 76 of 
Lecture Notes in Computational Science and Engineering, Springer, 2011, pp. 43-62. 


29 


